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The projected entangled-pair state (PEPS) ansatz can represent a thermal state in a strongly 
correlated system. We introduce a novel variational algorithm to optimize this tensor network 
whose essential ingredient is an auxiliary tree tensor network (TTN). Since full tensor environment 
is taken into account, with increasing bond dimension the PEPS-TTN ansatz provides the exact 
Gibbs state. Our presentation opens with a ID version for a matrix product state (MPS-TTN) and 
then generalizes to PEPS-TTN in 2D. Benchmark results in the quantum Ising model are presented. 
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I. INTRODUCTION 

Since their conception as the density matrix renormal¬ 
ization group (DMRG) |1] - an algorithm to optimize the 
matrix product state (MPS) ansatz in ID [2] - tensor net¬ 
works proved to be a competitive tool to study strongly 
correlated quantum systems. In the last decade, MPS 
was generalized to a 2D projected entangled pair state 
(PEPS) [3] and supplemented with an alternative multi¬ 
scale entanglement renormalization ansatz (MERA) [3]. 
These tensor networks avoid the notorious fermionic sign 
problem [5] and PEPS was applied to the t-J model of 
high-T c superconductivity providing the best results on 
the market 0. The networks - both MPS [2 [5] and 
PEPS roTTTTl - also made some major breakthroughs in 
search for topological order. 

Unlike the ground state, thermal states were explored 
mainly with MPS in ID [T31 ESj where they can be pre¬ 
pared by imaginary time evolution. One can follow simi¬ 
lar lines in 2D P1H51 the PEPS manifold is an efficient 
representation for Gibbs states na - but accurate evolu¬ 
tion is more demanding. Alternative direct contractions 
of the partition function were proposed [16] but - due to 
local tensor update - they are not warranted to become 
exact with increasing refinement parameter. In the fol¬ 
lowing we introduce a variational algorithm to optimize 
PEPS at finite temperature that both employs full tensor 
update and avoids direct imaginary time evolution. 


II. PURIFICATION OF THERMAL STATES 

We consider spins on an infinite lattice with a Hamil¬ 
tonian PL. Every spin has states s = 0,..., S — 1 and is 
accompanied by an ancilla with states a = 0, ...,S — 1. 
The enlarged “spin+ancilla” space is spanned by states 
n m | Sm,a, m ) where m numbers lattice sites. The Gibbs 
state of spins at an inverse temperature /3 is obtained 
from its purification |V’(/3)) in the enlarged space, 

( 1 ) 



FIG. 1. In A, the Trotter tensor To in ID. This spin operator 
with (red) spin indices depends on (black) bond indices con¬ 
necting it with similar operators at the nearest neighbor sites. 
In B, evolution operator U(/3) is a product of N time steps 
represented by horizontal rows of tensors. Each connecting 
line is a contracted index. In C, at every site the column 
of elementary To’s is compressed by an isometry W into one 
tensor T n with a bond dimension D < 2 N . In D, evolution 
operator U (/3) is approximated by a matrix product operator 
of the compressed tensors T n . 

At /3 = 0 we choose a product over lattice sites, 

S-i 

\m) = U ^2\ s m,s m ), ( 2 ) 

m s—0 

to initialize the imaginary time evolution, 

\m) = M0)> = U{P) \m)- (3) 

Here the Hermitian U(j3) = acts in the Hilbert 

space of spins. With the initial state ([2]) the trace in Eq. 
Q yields 

p(P) « U({3) U((3). (4) 

In the following the operator U(/3) will be represented 
by a projected entangled-pair operator (PEPO). Thanks 


p(P) °c Tr ancmas |^(/3))(V>(/3)|. 
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FIG. 2. The huge compression in Fig. HP is split into n 
steps. In A, as a first step k elementary Trotter tensors To 
are compressed by isometry Wi into tensor Ti with a bond 
dimension D < 2 . In B, the first step is followed by iterative 
compressions T m _i —> T m preserving the bond dimension D. 
In C, the final T n is the same as if the huge isometry W 
were a tree tensor network (TTN) with n layers of isometries: 
Wi ,..., W„. Here we show an example with k = 3 and n = 3. 


to the simple Eq. © equation [3]) translates between 
the PEPO and a PEPS for the purification \ip{f3)) in a 
trivial way: to obtain the PEPO it is enough to rebrand 
the PEPS’s ancilla indices as spin indices. 


III. QUANTUM ISING MODEL 



FIG. 3. In A, the partition function Z = (i!>(/3)\ip(/3)) is 
a product of transfer matrices t. In B, the same as in A 
but with one of the tensors T n removed. This diagram is 
the tensor environment E(n) for T n . In C, a lower-level en¬ 
vironment E(m — 1) for T m _i is obtained from E(m). In 
D, environment Ew(m) for isometry W m is obtained from 
E(m). In E, environment Ew{ 1) is obtained from E(l). In 
F, isometric environment is subject to singular value decom¬ 
position, Ew(m) = U AV , with D non-zero singular values 
A, a D 2 x D isometry U, and a D x D unitary V. The isome¬ 
try is updated as W rn = UV^ to maximize the figure of merit 
Z = Tr E w {m)W^ l . 


We proceed with the quantum Ising model: 

n = - ZmZ m' -h^Xm. (5) 

( 771 , 771 /) 771 

Here Z , A' are Pauli matrices and h is a transverse field. 
In ID, the model has a quantum critical point at h = 1 
that becomes a crossover at finite temperature. On a 2D 
square lattice, there is a ferromagnetic phase for small h 
and large /3. At zero temperature the quantum critical 
point is ho = 3.044, see Ref. |20j, and at h = 0 the 
Onsager’s critical point is /3o = — ln(v / 2 — l)/2 = 0.441. 
Lattice symmetries are not broken in any phase. 


IV. SUZUKI-TROTTER DECOMPOSITION 

We define gates 

U Z z{P)= 1] U x (/3) = l[e' hX -. (6) 

( 771 , 771 ') 771 

In the second-order Suzuki-Trotter decomposition a small 
time step is approximated by their product 

U(d0) « U x (dP/2)Uzz{d0)U x (d0/2). (7) 


In order to rearrange U(/3) as a tensor network, at each 
bond we make a singular value decomposition 

e ^rZ m z m , = ^ Zm b Zm ', b . ( 8 ) 

6 = 0,1 

Here b is a bond index and z mi b = V^b{Z m ) b with the 
singular values Ao = cosh /p and Ai = sinh Now we 
can write 



Here {b} is a set of all bond indices . The brackets 

enclose an elementary Trotter tensor To at site to. It is a 
spin operator depending on bond indices connecting its 
site with its nearest neighbors. Its ID version is shown 
in Fig. [TjA. 

The evolution is a product of N small time steps, 

U((3) = [umf- ( 10 ) 

For pedagogical reasons, we begin with a ID version of 
this product in Fig. [TjB, where each row l...iV is the 
elementary time step U (d/3) and each column is a site in 
the ID chain. We need an efficient algorithm to make 
this huge sum. 
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V. MPS-TTN ALGORITHM IN ID 

The general idea is to compress first the N elemen¬ 
tary tensors Tq at every site into a single tensor T n as in 
Fig. [Tp and then to contract the compressed T„’s hori¬ 
zontally as in Fig. HP- Here W is an isometry from the 
auxiliary Hilbert space spanned by 2 N possible values of 
N bond indices to its D-dimensional subspace. The bond 
dimension D is a refinement parameter: with increasing 
D results should become numerically exact. 

Furthermore, instead of making the huge compression 
W all at once - that is all but possible - we split it into 
n steps as in Figs. HMP- The final T n is the same 
as if W were a tree tensor network (TTN) [18] with n 
layers of isometries W\,...,W n , see the example in Fig. 
HP- Notice that the number of different isometries is only 
logarithmic in the number N = 2 n ~ 1 k of time steps d/3: 

N 0 

n = l + log 2 — ~ log AT = log-^r. (11) 
k dp 

Low temperatures or near-infinitesimal dj.3 required in 
precise applications can be achieved with a marginal log¬ 
arithmic cost. In the following we describe how this vari¬ 
ational TTN ansatz can be optimized in an efficient way. 

The optimization aims to maximize the partition func¬ 
tion 

Z [W] = e~ fiF = Tr e~ pH = Tr U{0)U{0) (12) 

with respect to the isometries Wi ,..., W„. This figure of 
merit has to be maximized in order to minimize the Gibbs 
free energy F. With each U(0) in Eq. Q represented by 
the diagram in Fig. [Tp), Z becomes the tensor network 
in Fig. [Sp¬ 
in order to maximize Z with respect to W m we need 
a gradient dZ/dW m . In principle, the gradient is an 
infinite sum of derivatives with respect to every W m in Z 
but, thanks to the symmetries, all these derivatives are 
the same and equal to a tensor environment %(m) of 
W m . The environment is the partition function in Fig. 
[3p, but with one tensor W m removed. 

The environment can be computed efficiently by an 
algorithm made of the procedures depicted in Figs[3p>- 
E. In Fig. [3j3 we show an environment E{n) for the 
compressed tensor T n . This environment is the parti¬ 
tion function in Fig[3j\, but with one tensor T n removed. 
It can be calculated with the standard transfer-matrix 
techniques applied to the transfer matrix t. Figure HP 
shows how to compute an environment E(m — 1) for ten¬ 
sor T m _i from a higher-level environment E{m) for T m . 
Repeating this procedure for m = n, ..., 2 we could obtain 
all Trotter tensors’ environments from E{n — 1) down to 
E(l) in one go. However, after each E(m) is obtained 
we pause to calculate an environment Ew(m) for the 
isometry W m by contracting the diagram in Fig. HP or 
its variant in Fig. [3p for m = 1. This environment is 
used immediately to optimize W m by the SVD technique 
depicted in Fig. HP- With the updated W rn the algo¬ 
rithm proceeds to calculation of E(m. — 1) and so forth 




FIG. 4. In A, energy per site of the quantum Ising chain 
at the critical field h = 1 in function of 0 for different bond 
dimensions D. For larger D the results remain accurate down 
to lower temperatures that are closer to the quantum critical 
point. In B, the correlation length in the tail of the ferromag¬ 
netic correlator (13). In both A and B, the number of isometry 
layers n = 11, the bond dimension D = 2 k — 2,..., 64, and 
the elementary time step in the Suzuki-Trotter decomposition 
d/3 = 0/(2 n ~ 1 k). 


all the way down to E( 1) and W\ when the down opti¬ 
mization sweep is finally completed. The end of the down 
sweep is the beginning of an up sweep whose first step 
is calculation of T) as in Fig. HP- With Ti we calculate 
E\v( 2) contracting the diagram in Figure HP and then 
immediately update W 2 before proceeding to calculation 
of T 2 . This procedures are repeated all the way up to 
W n and T n when the up sweep is completed. Before the 
next down sweep the environment E[n) in Fig. [3p is 
updated. The whole optimization procedure is repeated 
until convergence. 

In summary, the variational TTN is optimized by re¬ 
peated up- and down-sweeps. The up-sweep is a sequence 

To -A £]y( 1) —> T\ —y .... -A Ew{n) —y T n —> E{n) 
followed by the down sweep 

E(n) —y E\v (n/ — y E/n — 1) — y — — y E/ 1) — y E\y/ 1). 

Each environment E\y(m) is used immediately to update 
W m . Progress of the optimization is monitored with a set 
of figures of merit Z m = Tr Eyy/m/W^. When all Z m 
become the same within presumed numerical accuracy, 
then the isometries are accepted as converged. The nu¬ 
merical cost of the ID algorithm scales like (D 2 ) 3 , where 
D 2 is the bond dimension of the transfer matrix t. 


VI. BENCHMARK RESULTS IN ID 

We applied the algorithm to scan thermal states at the 
critical field h = 1 up to 0 ss 40, see Fig. [4] where we 
show the energy per site and the correlation length £ in 
the ferromagnetic correlator 

Cr = (Z X Z X+R ) - (Z X ){Z X+R ) 


(13) 
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FIG. 5. In A, the elementary Trotter tensor To on a square 
lattice. In B, compression of two tensors T m -i into one T m 
with four isometries W m . In C, tensors T n are contracted into 
a PEPS network for the evolution U(/3). In D, a contraction of 
two tensors T n makes a transfer tensor t. In E, contraction of 
transfer tensors is the partition function Z. In F, the partition 
function with one of the transfer tensors removed is a tensor 
environment Et. for the removed tensor t. We enclose the ends 
of the free bonds by a transparent reddish sphere. Each free 
bond has dimension D 2 . In G, when the environment/sphere 
E t is filled and contracted with one Trotter tensor T n we ob¬ 
tain tensor environment E(n) for T n . This environment is the 
partition function with one T n removed. In H, environment 
E(m) filled and contracted with two T m _Ts and three W m ’s 
makes environment Ewi'm) for isometry W m ■ In I, E(m — 1) 
is obtained from E{m). 


with an exponential tail Cr oc e~ R '^ for large R. We 
find that for increasing bond dimension D = 2 k results 
remain accurate down to lower temperatures, reaching 
closer to the quantum critical point. 

The number of TTN layers was fixed at n = 11 to 
achieve small enough d/3 for the results to be converged 
in d/3 all the way down to /3 ~ 40. Increasing n up to 20 
did not affect stability of the algorithm. The plots in Fig. 
[4] are scans in increasing j3 in the sense that isometries 


converged for one /3 were used as the initial isometries 
for the next /3 + 6/3. This recycling reduced the number 
of optimization sweeps necessary for each /3 to 10... 100. 
This number increased with /3 as the quantum critical 
point was approached. 

VII. PEPS-TTN ALGORITHM IN 2D 

The algorithm can be generalized to a 2D square lat¬ 
tice as summarized in Figure [5] The elementary Trotter 
tensor To now has four bond indices to be contracted with 
its four nearest neighbors, see Fig. HP- Accordingly, each 
compression of two Trotter tensors T m _ 1 into a higher 
level tensor T m is done with four isometries W rn , see Fig. 
[5j3. After all the isometric compressions are completed, 
the operator t/(/3) can be represented by the PEPO in 
Fig. [5p with tensors T n . A contraction of two TVs makes 
a transfer tensor t , see Fig. HP- The infinite network of 
transfer tensors in Fig. HP represents the partition func¬ 
tion Z. This 2D network replaces the ID chain of transfer 
matrices in Fig. HP- When one of the transfer tensors in 
Fig. §5 is removed from the partition function, then 
we obtain its tensor environment E t represented by the 
network in Fig. [5p. 

Unlike its analogue in ID, where the exact transfer- 
matrix techniques can be used, this environment requires 
more sophisticated technology. In this paper we use the 
symmetric version of the corner matrix renormalization 
(CMR) [19] whose description and discussion is delegated 
to appendix[A] Numerically it is the most expensive part 
of the algorithm with an additional refinement parameter 
of its own: the environmental bond dimension M. 

Once E t is converged, we can continue with the main 
loop of the algorithm. In Fig. HP E t is contracted with 
T n to yield an environment for T n that we call briefly 
E(n). With E{n) the down optimization sweep begins 
that proceeds as follows. From E{n) we obtain an en¬ 
vironment Ew{n) for the isometry W n as in Fig. [5p. 
Once calculated, this environment is used immediately 
to update W n as in Fig. HP- With the updated W n we 
proceed to calculate the environment E{n — 1) for T n _ 1 
as in Fig. [5j and then Ew{n — 1) to update W„_The 
same procedure repeats all the way down to W\. 

After all the isometries were optimized down to Wi, 
the upwards optimization sweep begins. It has n steps. 
In the m-th step two tensors TV-i and the environment 
E{m) - that was calculated before during the down sweep 
- are contracted to obtain Ew{m), see Fig. [5p. With this 
environment W m is updated immediately and then used 
to compress two T m -i into one T m as in Fig. [5p. This 
basic step is repeated all the way up to T n . 

Just as in ID, the above description can be briefly sum¬ 
marized as follows. The up-sweep is a sequence 

To —> Ew{ 1) —* Ti —> .... —> T n _ 1 —>• E\y{n ) —> T n 
completed with the CMR procedure 
T n ^ Et ^ E(n) . 
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FIG. 6. In A, spontaneous magnetization (Z) in the fer¬ 
romagnetic phase at the transverse field h = |/io- Different 
colors correspond to bond dimensions D = 2...8. Plots for 
D = 6, 7, 8 collapse demonstrating convergence for D > 6. In 
B, the ferromagnetic correlator Cr in Eq. (15 I for D = 6 at 
h = | ho and different (3 in the ferromagnetic phase. Its corre¬ 
lation length diverges closer to criticality requiring a diverging 
M. For M = 68 the longest length achieved is £ = 290 lattice 
sites. 




FIG. 7. In A, spontaneous magnetization (Z) near the crit¬ 
ical point. With increasing bond dimension D the plots be¬ 
come “critical” enough to make a fit (Z) oc (/3 — /3 C )' 3 with a 
critical-point j3 c = 0.5898 and a critical exponent /3' = 0.198. 
In B, the ferromagnetic correlation function ( |15| ) for D = 6 
at /3 C . In the intermediate range 1 < R < 30 the correlator 
is fitted well by the power law Cr oc R~ v with the critical 
exponent rj = 0.265. 


that is also the starting point for the down-sweep, 

E(n) —> Ew{n) —> E(n — 1 ) —> .... —> E{ 1 ) —> E\y{ 1 ) • 

Here, both in the up- and down-sweeps, each Ew{m) is 
used immediately to update W m . The whole up-down 
cycle is repeated until convergence. 

The numerical cost of the procedures in panels 
B,D,G,H,I of Fig. [H] scales formally like ( D 2 ) 4 , where 
D 2 is the bond dimension of the transfer tensor. The 
exponent is steeper than in ID, but a much smaller D is 
typically needed in 2D. For instance, unlike in ID, D can 
be finite at a critical point. However, the actual bottle¬ 
neck is the CMR in appendix [Aj whose cost scales like 
(D 2 ) 3 M 3 . In non-symmetric versions of CMR TS) this 
formal cost can be cut down to ( D 2 ) 2 M 3 . The parame¬ 
ter M is expected to diverge at a critical point together 
with a diverging correlation length but, even at critical¬ 
ity, local observables and correlations at increasingly long 
distance can be converged with increasing M. 


VIII. BENCHMARK RESULTS IN 2D 

We applied the PEPS-TTN algorithm to scan the fer¬ 
romagnetic phase at the transverse field 

h=^h 0 . (14) 

Results in Fig. [ 6 ] show convergence in the bond dimen¬ 
sion (for D > 6 ). For each D they are converged in the 
environmental bond dimension M . The plots do not ex¬ 
tend all the way down to the critical point, because the 
diverging correlation length would require a diverging en¬ 
vironmental bond dimension M to obtain a Gibbs state 
fully converged in M . Nevertheless, even with a limited 
M < 68 we could get close enough to the critical point 
to obtain a converged ferromagnetic correlation function 

C R = (Z X ^yZ X + R y) {Z X +R,y) (1^) 

with a long correlation length up to £ = 290 in its expo¬ 
nential tail Cr oc e~ R ^ for large R. 

The plots in Fig. [ 6 ] are fully converged in M, but 
they terminate before the critical point provoking a nat¬ 
ural question what, if anything, can be achieved closer 
to criticality. Hence we pushed our computations closer, 
accepting the fact that the ansatz cannot be fully con¬ 
verged there. Results are shown in Fig. [7J The spon¬ 
taneous magnetization for D = 6 and M = 35 allows a 
power law fit that gives the order parameter exponent 
P' = 0.198 (the textbook exponent P distinguished here 
by a prime from the inverse temperature) and the crit¬ 
ical point (3 C = 0.5898. The correlation function at P c 
has an exponential tail that is not converged in M - its 
range increases with M - but at an intermediate range 
1 < R < 30 the correlator has the correct form Cr oc R~ v 
with the critical exponent 77 = 0.265. 

To give an idea about practical effectiveness of the al¬ 
gorithm, it would take 1 day on a 4-core laptop computer 
to reproduce Fig. [ 6 ] and, since the convergence is slower 
near criticality, 2-3 more days for Fig. [TJ Our TTN had 
n = 6 layers of isometries with k = 5 in the bottom layer. 
The computation time was checked to be practically inde¬ 
pendent of n, as it is CMR that is the actual bottle-neck 
and not the isometry optimizations. For random initial 
conditions, the number of up- and down- optimization 
sweeps necessary to reach convergence was ~ 100. For 
a scan in P - when tensors converged for one P were re¬ 
cycled as initial tensors at a near p — SP - the number 
of sweeps was typically ~ 10. Both numbers increased 
towards the critical point. 


IX. COMPARISON WITH DIRECT 
IMAGINARY TIME EVOLUTION 


In our previous work [l4( 115] , we made efforts to ob¬ 
tain thermal states by direct imaginary time evolution 
from P = 0 to a finite one. In Fig. [8] we compare our 
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FIG. 8. Comparison between the present variational method 
in A and the direct imaginary time evolution in B. The TTN 
ansatz is exponentially more compact than MPS since the 
number of different isometries is only logarithmic in the num¬ 
ber of time steps. On top of this, in the present variational 
method the isometries in TTN are optimized by repeated up- 
and down-sweeps to provide only the most accurate state at 
the final /3, while in the time evolution all intermediate states 
between 0 and the final (3 need to be accurate compromising 
the accuracy of the final one. 


present method with the time evolution. The most strik¬ 
ing difference is the ansatz for the huge isometry W. 
In the time evolution we apply one elementary Trotter 
tensor To at every time step. Every application is fol¬ 
lowed by a renormalization of the PEPS tensor with a 
new isometry W m . This procedure implicitly assumes 
the (left-canonical) matrix product state (MPS) ansatz 
in Fig. [8j3. It has 0((3/d(3) different isometries, while 
in the corresponding TTN in Fig. the same number 
is only O [\og((3/d/3)\. This is a dramatic difference for 
large (3 or precise applications that require infinitesimal 
d(3. 

Another, possibly less apparent, difference is the op¬ 
timization procedure. In the time evolution it is a one 
way process. After the m-th time step we choose the 
optimal isometry W m to renormalize the bond indices of 



FIG. 9. Spontaneous magnetization (Z) in function of (3 ob¬ 
tained with different algorithms. The direct imaginary time 
evolution produces a plot with a discontinuity near the second 
order transition due to insufficient M. Adding a tiny sym¬ 
metry breaking bias 5 = 10“ 6 smooths the transition making 
the necessary M finite. Both evolution curves come from Ref. 
[Tol . They are compared with the PEPS-TTN curve from Fig. 
EK Here all the results for D = 6. 

the new PEPS tensor after this step. Once chosen, W m 
remains fixed and an error incurred with W m affects the 
whole following evolution. The errors accumulate with 
time. In the present variational approach we assume an 
entirely different strategy. Instead of scanning the whole 
range from 0 to /3 we target only the final (3. The isome¬ 
tries are optimized by repeated up- and down-sweeps to 
provide the best final state. Each isometry alone and 
the whole set of isometries in tune are serving the single 
goal to make the targeted state as accurate as possible. 
Unlike in the time evolution, they are not compromised 
to provide accurate thermal states at intermediate imag¬ 
inary times. Even if T n is the best PEPS tensor at the 
final f3, T n _ s does not need to be the best one at (3/s. 

The last property contributes to the main advantage 
of the variational PEPS-TTN over time evolution. The 
variational algorithm can probe a low temperature phase 
without any need to evolve from infinite temperature 
across a critical point. Not only there is no direct evo¬ 
lution through intermediate temperatures, but also the 
intermediate tensors T n _ s do not need to be thermal 
states at all. The problems with evolution across the crit¬ 
ical point are illustrated in Figure [9j where we compare 
three ferromagnetic magnetization curves obtained with 
three different algorithms. The direct evolution runs into 
trouble at the critical point, where exact evolution would 
require a divergent M, and the magnetization makes a 
discontinuous jump. This problem can be partly circum¬ 
vented by adding a tiny symmetry breaking bias to the 
Hamiltonian, 

A H = (16) 

m 

that smooths the transition making the necessary M fi¬ 
nite, but significantly alters the physics at the critical 
point. Nevertheless, the bias allows smooth evolution 
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from infinite temperature deep into the low temperature 
phase, where the effect of the tiny bias becomes negligi¬ 
ble. However, even the bias cannot prevent the evolution 
from accumulating errors with time. 

In retrospect, it may be tempting to combine the 
new variational strategy with the MPS ansatz instead of 
TTN. After all, the powerful methods developed for MPS 
[2] may be efficient enough to optimize even the huge 
number 0(/3/<i/3) of isometries. Unfortunately, these 
methods cannot be applied in 2D. Instead, the isometries 
in MPS have to be optimized by the procedures in Fig. 
[5] combined with CMR. Since an isometry in MPS maps 
from 2D to D dimensions, rather than from D 1 2 3 4 5 to D in 
TTN, the procedures are more efficient for MPS than for 
TTN. However, the actual bottleneck that limits D on an 
infinite lattice is the CMR whose cost depends only on 
the net D and not on the underlying ansatz. Thus with 
MPS one can achieve the same D as with TTN, but at 
the expense of an algorithm that is linear instead of log¬ 
arithmic in fi/dfi. Apart from this, the linear algorithm 
has a lot more variational parameters, hence in principle 
it may be more liable to getting trapped in local maxima 
of the figure of merit. However, this discussion does not 
quite exclude MPS in some applications like, e.g., a finite 


lattice. 


X. CONCLUSION 

The proposed PEPS-TTN method is the first 2D finite- 
temperature tensor-network algorithm that is both vari¬ 
ational and becomes numerically exact with increasing 
bond dimension. It avoids the demanding direct imagi¬ 
nary time evolution and employs full tensor environment 
in variational optimization. There is still a lot of room for 
improvement and development. For instance, the simple 
TTN could be made more powerful with some unitary 
disentanglers [4j. Internal symmetries, like Z 2 or U( 1), 
could help to use the bond dimension more efficiently. 
Finally, CMR may be upgraded to an algorithm making 
more efficient use of the environmental bond dimension. 
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FIG. 10. On the left, a planar version of the network in 
Fig. [5^ representing the partition function. Here each (black) 
bond represents two bond indices in Fig. [5j3. Its dimension is 
D 2 . This infinite contraction cannot be done exactly, hence 
it is approximated by the finite network on the right. The 
finite corner matrices C and top tensors T effectively repre¬ 
sent corresponding infinite sectors of the network on the left 
separated by the dashed blue lines. Their (red) environmen¬ 
tal bonds have dimension M. The environmental tensors C 
and T should be such that to the transfer tensor t in the cen¬ 
ter its effective environment on the right appears the same as 
its exact environment on the left as much as possible. They 
are obtained by iterating until convergence the corner matrix 
renormalization in Fig. 
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Appendix A: Corner matrix renormalization (CMR) 

An infinite tensor network, like the one in Fig. [5]E, 
cannot be contracted exactly. Fortunately, what we need 
in general is not this number, but an environment for 
a few tensors of interest. For instance, in Fig. [5jT we 
want an environment E t for the transfer tensor t in the 
center. The environment is a tensor that remains after 
removing the central tensor t from the infinite network. 
From the point of view of the central tensor, its infinite 
environment can be substituted with a finite effective en- 



FIG. 11. The corner C and top T are obtained by repeating 
until convergence a renormalization procedure. The proce¬ 
dure has four steps. In the first step, the tensors C, T, t are 
contracted to an enlarged corner C". In the second step, the 
symmetric MD 2 x MD 2 matrix C" is diagonalized and its M 
eigenvectors with the largest eigenvalues define an isometry 
Z. The diagonalization that scales like M 3 D 6 is the leading 
cost of this variant of corner matrix renormalization. In the 
third step, Z is used to renormalize/truncate the indices of 
C" back to the original dimension M giving a new (diagonal) 
corner C'. In the fourth step, the same Z renormalizes the 
contraction of T with t to a new T'. The four-step procedure 
is repeated until convergence of the M leading eigenvalues of 
C". 

B) 


FIG. 12. In A, a finite tensor network approximating the 
infinite environment Et in Fig. [5^. Notice that the corner 
matrix C can be singular-value-decomposed, C = UXV\ and 
then absorbed into the neighboring top tensors T. With X = 
s/XU^TUVX and Y = VXV^TV-\f\ (here the matrix products 
are understood in the environmental bond indices) we obtain 
an equivalent network in panel B. This network is a matrix 
product state (MPS) of four tensors X, Y, X, Y with a bond 
dimension M, hence in principle a finite M could suffice to 
represent the exact Et even at the critical point. 
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vironment, made of finite corner matrices C and top ten¬ 
sors T, that appears to the central tensor the same as the 
exact environment as much as possible. The environmen¬ 
tal tensors are contracted with each other by indices of 
dimension M. Increasing M makes the effective environ¬ 
ment more accurate and, for a finite correlation length, 
the effective environment is expected to converge to the 
exact one at a finite M. In the Ising model, tensor t is 
symmetric under permutation of its indices and, conse¬ 
quently, C is a symmetric matrix and T is symmetric in 
its environmental indices. 

Finite tensors C and T represent infinite sectors of the 
network on the left of Fig. 10 The tensors are con- 
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verged by iterating the corner matrix renormalization in 
Fig. [TT} In every renormalization step, the corner ma¬ 
trix is enlarged with one t and two T’s. This operation 
represents the top-left corner sector in Fig. [T0| absorbing 
one more layer of tensors t. Once the environment is con¬ 
verged, it can be used to calculate either an observable 
or the environment E t , see Fig. |12| 

The above CMR procedure requires M that diverges 
at a critical point. This is clearly demonstrated in the 
appendix of Ref. |15] at the Onsager transition in the 2D 
classical Ising model, where a finite M results in a finite 
correlation length oc M 193 . The tail of the correlation 
function cannot become strictly algebraic for any finite 
M but, as illustrated in section VII and Ref. [15] , with 
increasing M not only local observables but also corre¬ 
lations at increasingly long distance become converged 
in M. At a critical point the effective environment is 
not exact, but provides an approximation whose quality 
improves with M in a systematic way. 

In an attempt to go beyond this bottom line, one could 


argue that the divergent M required at criticality is not 
an inherent property of the finite effective environment, 
but an artifact of the specific CMR method used to ob¬ 
tain this environment. Indeed, the exact tensor E t in Fig. 
[5^ can be interpreted as a quantum state on the four free 
bonds, each free bond with D 2 auxiliary states. Such a 
finite state can be represented by a four-site MPS with 
a finite bond dimension M. In Fig. [12) 3 we show that 
this MPS is equivalent to a finite effective environment 
with the same M. This completes the argument that for 
a finite D the exact E t can be represented with a finite 
M. 

An exact local E t is all we need to optimize the isome¬ 
tries W m of the thermal PEPS, even though the finite M 
of its exact environmental tensors may prohibit accurate 
calculation of the power-law tails of critical correlations. 
The present CMR method attempts to be universal - ac¬ 
curate for both local and non-local observables - while 
in order to optimize the PEPS tensor all we need is a 
method targeting the local E t only. 



